!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!
!
C
C   /* Deck hucdrv */
      SUBROUTINE HUCDRV(CMO,WORK,LWORK)
C
C     K.Ruud, July-96
C
#include "implicit.h"
#include "maxorb.h"
C
      DIMENSION WORK(LWORK), CMO(NCMOT)
#include "infdim.h"
#include "inforb.h"
#include "huckel.h"
C
      CALL QENTER('HUCDRV')
      NROW = NBAST + NHUCBA
      NELMNT = NROW*(NROW + 1)/2
      NHUCMX = 0
      DO 10 ISYM = 1, NSYM
         NHUCMX = MAX(NHUCMX,NHUCAO(ISYM))
 10   CONTINUE
      NNHUCMX = NHUCMX*(NHUCMX + 1)/2
      N2HUCMX = NHUCMX*NHUCMX
      N2HUCBA = NHUCBA*NHUCBA
      NTRDIM = MAX(NHUCMX*NBASMA,NBASMA*(NBASMA + 1)/2)
C
      IF (NHUCBA .LE. 0) THEN
         CALL QUIT('No Huckel basis functions.')
      END IF
C
      KFREE = 1
      LFREE = LWORK
      CALL MEMGET('REAL',KTOTOV,NELMNT,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KHUCMT,NNHUCMX,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KHUCOV,N2HUCBA,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KEIGVC,N2HUCMX,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KTRANM,NTRDIM,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KEIGVL,NHUCBA,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KEIGVS,NHUCBA,WORK,KFREE,LFREE)
      CALL MEMGET('REAL',KOVERL,NBASMA*NBASMA,WORK,KFREE,LFREE)
      CALL HUCDRV_1(WORK(KFREE),LFREE,WORK(KTOTOV),NELMNT,NROW,
     &            WORK(KHUCMT),WORK(KHUCOV),NHUCMX,WORK(KEIGVC),
     &            WORK(KTRANM),WORK(KEIGVL),WORK(KEIGVS),WORK(KOVERL),
     &            CMO)
      CALL MEMREL('HUCDRV',WORK,1,1,KFREE,LFREE)
      CALL QEXIT('HUCDRV')
      RETURN
      END
C   /* Deck HUCDRV_1 */
      SUBROUTINE HUCDRV_1(WORK,LWORK,TOTOVL,NELMNT,NROW,
     &                  HUCMAT,HUCOVL,NHUCMX,
     &                  EIGVEC,TRANMT,EIGVAL,EIGVLS,OVERLP,CMO)
C
C     Driver routine for solving the extended Huckel problem, and 
C     projecting the eigenvectors onto the complete basis set
C     
C     K.Ruud, July-96
C
#include "implicit.h"
#include "dummy.h"
#include "maxorb.h"
      LOGICAL ANTSYM, FNDLB2
      CHARACTER*8 RTNLBL(2)
      DIMENSION WORK(LWORK)
      DIMENSION TOTOVL(NELMNT), HUCMAT(*), HUCOVL(*), TRANMT(*),
     &          EIGVEC(NHUCMX*NHUCMX), CMO(NCMOT),
     &          EIGVAL(NHUCBA), EIGVLS(NHUCBA), OVERLP(NBASMA*NBASMA)
C
#include "priunit.h"
#include "inftap.h"
#include "infdim.h"
#include "inforb.h"
#include "infpri.h"
#include "scbrhf.h"
#include "huckel.h"
C
      IF (P6FLAG(6) .OR. IPRI6.GE.0 .OR. IPRRHF.GE.0 .OR. EWMO)
     &   WRITE (LUPRI,'(/A/2(/5X,A,L10)/5X,A,8I5)')
     &      ' *** Output from Huckel module :',
     &      'Using EWMO model: ',EWMO,
     &      'Using EHT  model: ',.NOT.EWMO,
     &      'Number of Huckel orbitals each symmetry:',
     &      (NHUCAO(I),I=1,NSYM)
      IF (EWMO) WRITE (LUPRI,'(/A/A/A)')
     &  ' EWMO - Energy Weighted Maximum Overlap - is a Huckel type'//
     &  ' method,',
     &  '        which normally is better than Extended Huckel Theory.',
     &  ' Reference: Linderberg and Ohrn, Propagators in Quantum'//
     &  ' Chemistry (Wiley, 1973)'
C
      NNHUCBA = NHUCBA * (NHUCBA + 1) / 2
C
C     Read AO overlap matrix from AOPROPER file
C
      CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ','UNFORMATTED',IDUMMY,
     &            .FALSE.)
      IF (EWMO) THEN
C
C     ... read EWMO overlap matrix:
C
         REWIND LUPROP
         IF (FNDLB2('HUCKOVLP',RTNLBL,LUPROP)) THEN
            CALL READT(LUPROP,NELMNT,TOTOVL)
            IF (P6FLAG(6) .OR. IPRRHF .GE. 10) THEN
               CALL HEADER(
     &         'HUCDRV_1: HUCKOVLP matrix read from AOPROPER',-1)
               CALL OUTPAK(TOTOVL,NROW,1,LUPRI)
            END IF
         ELSE
            WRITE (LUPRI,'(//3A)') ' HUCDRV_1 error: "HUCKOVLP"'//
     *         ' matrix not found on AOPROPER.'
            CALL QUIT('HUCDRV_1 error: HUCKOVLP not found on AOPROPER')
         END IF
C
C     ... extract EWMO overlap matrix from TOTOVL:
C
         CALL DZERO(HUCOVL,NNHUCBA)
         ISTRT  = 0
         NHUCT  = 0
         JPOS   = 1
C
         DO ISYM = 1, NSYM
            NHUCAOI = NHUCAO(ISYM)
         IF (NHUCAOI .GT. 0) THEN
            NBASI  = NBAS(ISYM)
            IBASTO = IBAS(ISYM) + NHUCT
            IHUCI  = IBASTO + NBASI
            ISTRT  = (IHUCI*(IHUCI + 1))/2
C
C           Copy the Huckel overlap matrix into a lower triangular matrix,
C           where the elements are sorted in consecutive order
C
            IPOS = JPOS
            DO I = 1, NHUCAOI
               CALL DCOPY(I,TOTOVL(ISTRT + IHUCI + 1),1,HUCOVL(IPOS),1)
               IPOS = IPOS + I
               ISTRT = ISTRT + (IHUCI + I)
            END DO   ! I = 1,NHUCAOI
            IF (P6FLAG(6) .OR. IPRRHF.GE.6) THEN
               WRITE (LUPRI,'(/A,I3,/A)')
     &            ' Extracted EWMO overlap matrix of symmetry :',ISYM,
     &            ' ------------------------------------------------'
               CALL OUTPAK(HUCOVL(JPOS),NHUCAOI,1,LUPRI)
            END IF
            JPOS = IPOS
            NHUCT = NHUCT + NHUCAOI
         END IF
         END DO   ! ISYM = 1,NSYM
      END IF
C
C     Read Huckel matrix, both for EWMO and EHT
C
      REWIND LUPROP
      IF (FNDLB2('HUCKEL  ',RTNLBL,LUPROP)) THEN
         CALL READT(LUPROP,NELMNT,TOTOVL)
         IF (P6FLAG(6)) THEN
            CALL HEADER(
     &      'HUCDRV_1: Huckel matrix read from AOPROPER',-1)
            CALL OUTPAK(TOTOVL,NROW,1,LUPRI)
         END IF
      ELSE
         WRITE (LUPRI,'(//3A)') ' HUCDRV_1 error: "HUCKEL  "'//
     &      ' matrix not found on AOPROPER.'
         CALL QUIT('HUCDRV_1 error:'//
     &      ' HUCKEL overlap not found on AOPROPER')
      END IF
      CALL GPCLOSE(LUPROP,'KEEP')
C
C     Save initial CMO for HUCORT
C
      KCMO1 = 1
      KWRK  = KCMO1 + NCMOT
      LWRK  = LWORK - (KWRK-1)
      CALL DCOPY(NCMOT,CMO,1,WORK(KCMO1),1)
C
C
C     Extract Huckel matrix, do one symmetry at a time
C
      ISTRT  = 0
      NHUCT  = 0
      JPOS   = 1
C
      DO 100 ISYM = 1, NSYM
         NHUCAOI = NHUCAO(ISYM)
       IF (NHUCAOI .EQ. 0) GO TO 100
         NBASI  = NBAS(ISYM)
         IBASTO = IBAS(ISYM) + NHUCT
         IHUCI  = IBASTO + NBASI
         NORBI  = NORB(ISYM)
         IF (NORBI .EQ. 0) THEN
            WRITE (LUPRI,'(/A,I3/I5,A/A)')
     &         'HUCDRV warning in symmetry',ISYM,
     &         NHUCAOI,' Huckel orbitals, but no orbitals!',
     &         'Thus no Huckel projection in this symmetry!'
            GO TO 100
c           CALL QUIT('HUCDRV error, NHUCAOI.ne.0 .and. NORBI.eq.0')
         END IF
         ISTRT  = (IHUCI*(IHUCI + 1))/2
         ICSYM  = ICMO(ISYM) + 1
C
C     Extract ordinary overlap matrix
C
         IOVLST = (IBASTO*(IBASTO + 3))/2 + 1
         IPOS = 1
         DO 115 I = 1, NBASI
            CALL DCOPY(I,TOTOVL(IOVLST),1,TRANMT(IPOS),1)
            IPOS = IPOS + I
            IOVLST = IOVLST + (IBASTO + I)
 115     CONTINUE
         CALL DSPTSI(NBASI,TRANMT,OVERLP)
C
C     Copy the Huckel part into a lower triangular matrix, where the elements
C     are sorted in consecutive order
C
C     Also extract the transformation matrix from HUCKEL basis to full basis
C
         IPOS = 1
         DO 110 I = 1, NHUCAOI
            CALL DCOPY(I,TOTOVL(ISTRT + IHUCI + 1),1,
     &                 HUCMAT(IPOS),1)
            CALL DCOPY(NBASI,TOTOVL(ISTRT + IBASTO + 1),1,
     &                 TRANMT((I - 1)*NBASI + 1),1)
            IPOS = IPOS + I
            ISTRT = ISTRT + (IHUCI + I)
 110     CONTINUE 
         IF (P6FLAG(6) .OR. IPRRHF.GE.8) THEN
            WRITE (LUPRI,'(/A,I3,/A)')
     &           ' Extracted Huckel matrix of symmetry :',ISYM,
     &           ' ------------------------------------------------'
            CALL OUTPAK(HUCMAT,NHUCAOI,1,LUPRI)
         END IF
         IF (P6FLAG(6) .OR. IPRRHF.GE.15) THEN
            WRITE (LUPRI,'(/A,I3,/A)')
     &           ' Extracted transformation matrix of symmetry :',ISYM,
     &           ' ------------------------------------------------'
            CALL OUTPUT(TRANMT,1,NBASI,1,NHUCAOI,NBASI,NHUCAOI,1,LUPRI)
         END IF
C
C       ??? EWMO or EHT ?
C
        IF (EWMO) THEN
C        ... call EISPACK routine RSG
         N2HUCAOI = NHUCAOI*NHUCAOI
         IPOS = 1 + (NHUCT*(NHUCT + 1))/2 + NHUCT
         KHMAT = KWRK
         KSMAT = KHMAT+ N2HUCAOI
         KFV1  = KSMAT+ N2HUCAOI
         KFV2  = KFV1 + NHUCAOI
         KLAST = KFV2 + NHUCAOI
         IF (KLAST .GT. LWORK) CALL STOPIT('HUCDRV','RSG',KLAST,LWORK)
         CALL DSPTSI(NHUCAOI,HUCMAT,WORK(KHMAT))
         CALL DSPTSI(NHUCAOI,HUCOVL(JPOS),WORK(KSMAT))
         IF (P6FLAG(6) .OR. IPRRHF.GE.8) THEN
            WRITE (LUPRI,*) 'EWMO AO Hamiltonian matrix W'//
     &                      ' corresponding to separate atoms'
            CALL OUTPUT(WORK(KHMAT),1,NHUCAOI,1,NHUCAOI,
     &                  NHUCAOI,NHUCAOI,1,LUPRI)
         END IF
         CALL DGEMM('N','N',NHUCAOI,NHUCAOI,NHUCAOI,1.D0,
     &              WORK(KHMAT),NHUCAOI, WORK(KSMAT),NHUCAOI,0.D0,
     &              EIGVEC,NHUCAOI)
         CALL DGEMM('N','N',NHUCAOI,NHUCAOI,NHUCAOI,1.D0,
     &              WORK(KSMAT),NHUCAOI, EIGVEC,NHUCAOI,0.D0,
     &              WORK(KHMAT),NHUCAOI)
C        HMAT = SWS , SMAT = S
         IF (P6FLAG(6) .OR. IPRRHF.GE.4) THEN
            WRITE (LUPRI,'(/A)') 'EWMO Hamiltonian matrix SWS'
            CALL OUTPUT(WORK(KHMAT),1,NHUCAOI,1,NHUCAOI,
     &                  NHUCAOI,NHUCAOI,1,LUPRI)
            WRITE (LUPRI,'(/A)') 'EWMO overlap matrix S'
            CALL OUTPUT(WORK(KSMAT),1,NHUCAOI,1,NHUCAOI,
     &                  NHUCAOI,NHUCAOI,1,LUPRI)
         END IF
         JPOS = JPOS + ( NHUCAOI * (NHUCAOI+1) ) / 2
         CALL RSG(NHUCAOI,NHUCAOI,WORK(KHMAT),WORK(KSMAT),
     &            EIGVAL(NHUCT+1),1,EIGVEC,WORK(KFV1),WORK(KFV2),IERR)
         IF ( IERR.NE.0 ) THEN
            WRITE(LUPRI,'(/A,I10)')
     *      ' HUCDRV EWMO eigenvalue problem not converged, code: ',IERR
            CALL QUIT('HUCDRV: EWMO eigenvalue equation not converged')
         END IF
C           Save orbital symmetry in EIGVLS
            DO I = 1, NHUCAOI
               EIGVLS(NHUCT + I) = ISYM
            END DO
        ELSE ! not EWMO
C
C       EHT: Diagonalize Huckel matrix (Note: KLAST is not a misprint)
C
         KSCRA1 = KWRK
         KSCRA2 = KSCRA1 + NHUCAOI
         KLAST  = KSCRA1 + NBASI*NHUCAOI
         IF (KLAST .GT. LWORK) CALL STOPIT('HUCDRV','JACO',KLAST,LWORK)
         CALL DUNIT(EIGVEC,NHUCAOI)
         CALL JACO(HUCMAT,EIGVEC,NHUCAOI,NHUCAOI,NHUCAOI,WORK(KSCRA1),
     &             WORK(KSCRA2))
C
C     Sort the eigenvalues and eigenvectors
C     Save orbital symmetry in EIGVLS
C
         II = 0
         DISYM = ISYM
         DO 120 I = 1, NHUCAOI
            II = II + I
            EIGVAL(NHUCT + I) = HUCMAT(II)
            EIGVLS(NHUCT + I) = DISYM
 120     CONTINUE 
         CALL ORDER(EIGVEC,EIGVAL(NHUCT + 1),NHUCAOI,NHUCAOI)
        END IF
C
C
C     We write the eigenvalues if high enough print level
C     and also the vectors if even higher print level
C
         IF (P6FLAG(6) .OR. IPRI6.GE.0 .OR. IPRRHF.GE.0) THEN
            IF (EWMO) THEN
               WRITE (LUPRI,'(/A,I3)') ' Huckel EWMO eigenvalues for'//
     &           ' symmetry :',ISYM
            ELSE
               WRITE (LUPRI,'(/A,I3)') ' Huckel EHT eigenvalues for'//
     &           ' symmetry :',ISYM
            END IF
            WRITE (LUPRI,'(5X,5F15.6)') (EIGVAL(NHUCT+I), I=1, NHUCAOI)
         END IF
         IF (P6FLAG(6) .OR. IPRRHF.GE.2) THEN
            IF (EWMO) THEN
               WRITE (LUPRI,'(/A,I3)') ' Huckel EWMO eigenvectors for'//
     &           ' symmetry :',ISYM
            ELSE
               WRITE (LUPRI,'(/A,I3)') ' Huckel EHT eigenvectors for'//
     &           ' symmetry :',ISYM
            END IF
            CALL OUTPUT(EIGVEC,1,NHUCAOI,1,NHUCAOI,
     &                  NHUCAOI,NHUCAOI,1,LUPRI)
         END IF
         NHUCT = NHUCT + NHUCAOI
C         
C     Project eigenvectors onto complete basis
C
C     We do this as (S_{NN})^{-1}*(S_{NO}*C_{O})
C     N and O refers to new (large) and old (small) basis respectively
C     For HUCKEL basis to HUCKEL basis this gives the identity matrix

C     Jan 2007 HJAaJ: No! We want to project on the MO basis
C     from DELMO, otherwise this will reintroduce some parts of the
C     deleted atomic orbitals again, leading to "mysterious" too
C     big MO coefficients.
C     Revised algorithm: (CMO_{Nn})^{T}*(S_{NO}*CMO_{OO})
C     where N,n is new basis NBASI,NORBI; O is old basis (Huckel basis);
C     and CMO_(OO) contains the Huckel molecular orbitals in Huckel basis.
C
C     ( Recall that CMO^{T} = CMO^{-1} * S^{-1} )

         KSCRA1 = KWRK
         KSCRA2 = KSCRA1 + NBASI*NHUCAOI
         KLAST  = KSCRA2 + NBASI*NHUCAOI
         IF (KLAST.GT.LWORK) CALL STOPIT('HUCDRV','DGEMM ',KLAST,LWORK)
C        1) Do (S_{NO}*CMO_{O}) part: 
         CALL DGEMM('N','N',NBASI,NHUCAOI,NHUCAOI,1.D0,
     &              TRANMT,NBASI, EIGVEC,NHUCAOI, 0.D0,
     &              WORK(KSCRA1),NBASI)
C        2) Check for orbital types in Huckel basis not in regular basis
C           (e.g. a 3d on Aluminium, not included in basis sets w/o pol)
         MHUCAOI = 0
         DO IHUCAOI = 1, NHUCAOI
            JSCRA1 = KSCRA1 + (IHUCAOI-1)*NBASI
            TNORM = DNRM2(NBASI,WORK(JSCRA1),1)
            IF (TNORM .GT. 1.0D-8) THEN
               JSCRA2 = KSCRA2 + MHUCAOI*NBASI
               MHUCAOI = MHUCAOI + 1
               CALL DCOPY(NBASI,WORK(JSCRA1),1,WORK(JSCRA2),1)
            END IF
         END DO
         IF (MHUCAOI .GT. NORBI) MHUCAOI = NORBI
C        ... hjaaj sep 07: decrease if huckel basis bigger than user basis,
C            e.g. if user is calculating with S-atom without d-orbitals.
         IF (NHUCAOI .NE. MHUCAOI) THEN
            WRITE(LUPRI,'(/A,I2,A,I0,A,I0)')
     &      'HUCDRV: reduced number of huckel orbitals in sym',ISYM,
     &      ' from ',NHUCAOI,' to ',MHUCAOI
            NHUCAOI = MHUCAOI
         END IF
C        3) Do (CMO_{Nn})^{T}* (work(kscra2)) part
         CALL DGEMM('T','N',NORBI,NHUCAOI,NBASI,1.D0,
     &              CMO(ICSYM),NBASI, WORK(KSCRA2),NBASI,0.D0,
     &              WORK(KSCRA1),NORBI)
         CALL DGEMM('N','N',NBASI,NHUCAOI,NORBI,1.D0,
     &              WORK(KCMO1-1+ICSYM),NBASI, WORK(KSCRA1),NORBI,0.D0,
     &              CMO(ICSYM),NBASI)
 100  CONTINUE   ! DO 100 ISYM = 1, NSYM
C
C     Determine Hartree-Fock occupation if requested
C
      IF (AUTOCC) THEN
         CALL ORDER (EIGVLS,EIGVAL,NHUCT,1)
         CALL IZERO(NRHF,8)
         MOCC = NRHFEL/2
         DO 98 IK = 1, MOCC
            ISYM = NINT(EIGVLS(IK))
            NRHF(ISYM) = NRHF(ISYM) + 1
 98      CONTINUE
         CALL ICOPY(8,NRHF,1,NOCC,1)
         IF (2*MOCC .NE. NRHFEL) THEN
            IOPRHF = NINT(EIGVLS(MOCC + 1))
            LSYM   = IOPRHF
            NOCC(IOPRHF) = NOCC(IOPRHF) + 1
         END IF
      END IF
C
C     We do our own Gram-Schmidt Orthogonalization of the trial vectors
C
      IF (P6FLAG(6)) THEN
         WRITE (LUPRI,'(/A)')
     &   ' Huckel guess MOs before Gram-Schmidt orthogonalization:'
         CALL PRORB(CMO,.FALSE.,LUPRI)
      END IF
      CALL HUCORT(CMO,WORK(KCMO1),TOTOVL,WORK(KWRK),LWRK)
      IF (P6FLAG(6)) THEN
         WRITE (LUPRI,'(/A)')
     &   ' Huckel guess MOs after Gram-Schmidt orthogonalization:'
         CALL PRORB(CMO,.FALSE.,LUPRI)
      END IF
C
C     We are done
C
      RETURN
      END
C   /* Deck hucfun */
      SUBROUTINE HUCFUN(NQ,IQCORE,IQM,ISHL1,HUCTER)
C
C panor/johhe 2005:
C     Added Hueckel start for elements up to radon (lanthanoids excluded) 
C     both for all-electron and effective core calculations. Ionizations
C     energies are taken from: Joseph B. Mann, "Atomic Structure Calculations 
C     I. Hartree-Fock Energy Results for the Elements H through Lr", 
C     Report LA-3690 (Los Alamos National Laboratory, 1967).
C
C hjaaj mar 2006:
C     Added 2p for Li, Be; 3p for Na, Mg; 4p for K, Ca
C     with energy parameters from A. Herman, Modelling Simul. Mater.
C     Sci. Eng. 12 (2004) 21. TODO: add 4p for first row transition metals etc.
C
C
C     Input:
C     NQ     = atomic number
C     IQCORE = number of core electrons described with ECP
C     IQM    = l_quantum_number + 1
C     ISHL1  = orbital no. for this IQM
C
C     Output:
C     HUCTER = orbital energy
C
#include "implicit.h"
#include "priunit.h"
#include "maxorb.h"
#include "huckel.h"
C
C
C     Ionization energies for shells:
C     1s,2s,3s,4s,5s,6s, 2p,3p,4p,5p,6p, 3d,4d,5d, 4f
C
      DIMENSION ENERGY(15,86), IOFF(4)
C
C     Ionization energies for 7. period shells:
C     1s,2s,3s,4s,5s,6s,7s,2p,3p,4p,5p,6p,7p,3d,4d,5d,6d,4f,5f
C
      DIMENSION ENERGY7(19,87:103), IOFF7(4)
C
      SAVE ENERGY, IOFF, ENERGY7, IOFF7
C
      DATA (IOFF(I),I=1,4)  /0,6,11,14/
      DATA (IOFF7(I),I=1,4) /0,7,13,17/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=1,10) /
C    H    Z= 1
     & 0.500000D+00,0.125000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.125000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   He    Z= 2
     & 0.918000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Li    Z= 3
     & 2.477700D+00,0.196300D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.130000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Be    Z= 4
     & 4.732700D+00,0.309300D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.206000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    B    Z= 5
     & 7.695300D+00,4.947000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,3.099000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    C    Z= 6
     & 11.33840D+00,7.121000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,4.069000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    N    Z= 7
     & 15.66640D+00,9.637000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,5.087000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    O    Z= 8
     & 20.68090D+00,1.250400D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.616200D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    F    Z= 9
     & 26.38270D+00,1.572500D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.730000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ne    Z=10
     & 32.77240D+00,1.930400D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.850400D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00/

C hjaaj Jun 2006: inserted Hydrogen 3d energy for Na-Ca
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=11,20) /
C   Na    Z=11
     & 40.47850D+00,2.797000D+00,0.182100D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,1.518100D+00,0.114000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Mg    Z=12
     & 49.03170D+00,3.767700D+00,0.253100D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,2.282200D+00,0.157000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Al    Z=13
     & 58.50090D+00,4.910700D+00,3.934000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,3.218300D+00,2.100000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Si    Z=14
     & 6.881920D+01,6.163000D+00,5.436000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,4.262400D+00,2.787000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    P    Z=15
     & 7.998890D+01,7.528700D+00,7.064000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,5.418400D+00,3.506000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    S    Z=16
     & 92.01040D+00,9.009900D+00,8.827000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,6.688000D+00,4.263000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Cl    Z=17
     & 1.048839D+02,1.060740D+01,1.072900D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,8.072200D+00,5.064000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ar    Z=18
     & 1.186098D+02,1.232210D+01,1.277300D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,9.571400D+00,5.910000D-01,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    K    Z=19
     & 1.335325D+02,1.448990D+01,1.748800D+00,1.475000D-01,0.000000D+00,
     & 0.000000D+00,1.151930D+01,9.544000D-01,0.103000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ca    Z=20
     & 1.493631D+02,1.682270D+01,2.245400D+00,1.955000D-01,0.000000D+00,
     & 0.000000D+00,1.362930D+01,1.340700D+00,0.138000D+00,0.000000D+00,
     & 0.000000D+00,0.055556D+00,0.000000D+00,0.000000D+00,0.000000D+00/

      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=21,30) /
C   Sc    Z=21
     & 1.658992D+02,1.908060D+01,2.567300D+00,2.101000D-01,0.000000D+00,
     & 0.000000D+00,1.566820D+01,1.574500D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,3.437000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ti    Z=22
     & 1.832827D+02,2.143490D+01,2.882000D+00,2.219000D-01,0.000000D+00,
     & 0.000000D+00,1.780300D+01,1.802900D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,4.058000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C    V    Z=23
     & 2.015220D+02,2.389690D+01,3.198800D+00,2.324000D-01,0.000000D+00,
     & 0.000000D+00,2.004440D+01,2.033500D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,4.612000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Cr    Z=24
     & 2.204413D+02,2.626840D+01,3.332400D+00,2.075000D-01,0.000000D+00,
     & 0.000000D+00,2.219840D+01,2.096200D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,3.225000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Mn    Z=25
     & 2.405799D+02,2.915980D+01,3.551000D+00,2.513000D-01,0.000000D+00,
     & 0.000000D+00,2.486240D+01,2.511100D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,5.612000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Fe    Z=26
     & 2.614001D+02,3.196460D+01,4.189200D+00,2.601000D-01,0.000000D+00,
     & 0.000000D+00,2.744260D+01,2.760500D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,6.079000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Co    Z=27
     & 2.830824D+02,3.488660D+01,4.536600D+00,2.686000D-01,0.000000D+00,
     & 0.000000D+00,3.013830D+01,3.017700D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,6.531000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ni    Z=28
     & 3.056269D+02,3.792650D+01,4.893600D+00,2.768000D-01,0.000000D+00,
     & 0.000000D+00,3.295030D+01,3.283000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,6.971000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Cu    Z=29
     & 3.287927D+02,4.081890D+01,5.011900D+00,2.385000D-01,0.000000D+00,
     & 0.000000D+00,3.561790D+01,3.324800D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,4.912000D-01,0.000000D+00,0.000000D+00,0.000000D+00,
C   Zn    Z=30
     & 3.533042D+02,4.436160D+01,5.637800D+00,2.925000D-01,0.000000D+00,
     & 0.000000D+00,3.892480D+01,3.839300D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,7.825000D-01,0.000000D+00,0.000000D+00,0.000000D+00/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=31,40) /
C   Ga    Z=31
     & 3.788181D+02,4.816830D+01,6.394600D+00,4.246000D-01,0.000000D+00,
     & 0.000000D+00,4.249390D+01,4.482300D+00,2.085000D-01,0.000000D+00,
     & 0.000000D+00,1.193400D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Ge    Z=32
     & 4.052501D+02,5.215600D+01,7.196800D+00,5.571000D-01,0.000000D+00,
     & 0.000000D+00,4.624190D+01,5.167300D+00,2.693000D-01,0.000000D+00,
     & 0.000000D+00,1.640600D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   As    Z=33
     & 4.326014D+02,5.632510D+01,8.044600D+00,6.951000D-01,0.000000D+00,
     & 0.000000D+00,5.016900D+01,5.895500D+00,3.301000D-01,0.000000D+00,
     & 0.000000D+00,2.127600D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Se    Z=34
     & 4.608717D+02,6.067350D+01,8.936700D+00,8.402000D-01,0.000000D+00,
     & 0.000000D+00,5.427350D+01,6.666100D+00,3.925000D-01,0.000000D+00,
     & 0.000000D+00,2.654200D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Br    Z=35
     & 4.900596D+02,6.519980D+01,9.871900D+00,9.927000D-01,0.000000D+00,
     & 0.000000D+00,5.855410D+01,7.478200D+00,4.571000D-01,0.000000D+00,
     & 0.000000D+00,3.220200D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Kr    Z=36
     & 5.201645D+02,6.990290D+01,1.084940D+01,1.152900D+00,0.000000D+00,
     & 0.000000D+00,6.300960D+01,8.331500D+00,5.242000D-01,0.000000D+00,
     & 0.000000D+00,3.825200D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Rb    Z=37
     & 5.514565D+02,7.504900D+01,1.213320D+01,1.523500D+00,1.377000D-01,
     & 0.000000D+00,6.790600D+01,9.487700D+00,8.101000D-01,0.000000D+00,
     & 0.000000D+00,4.732300D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Sr    Z=38
     & 5.836870D+02,8.039050D+01,1.347500D+01,1.896800D+00,1.785000D-01,
     & 0.000000D+00,7.299580D+01,1.069990D+01,1.098200D+00,0.000000D+00,
     & 0.000000D+00,5.694400D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C    Y    Z=39
     & 6.167480D+02,8.581060D+01,1.475890D+01,2.168900D+00,1.961000D-01,
     & 0.000000D+00,7.816420D+01,1.185420D+01,1.301200D+00,0.000000D+00,
     & 0.000000D+00,6.599500D+00,2.498000D-01,0.000000D+00,0.000000D+00,
C   Zr    Z=40
     & 6.507100D+02,9.138430D+01,1.606170D+01,2.424700D+00,2.086000D-01,
     & 0.000000D+00,8.348520D+01,1.302650D+01,1.492500D+00,0.000000D+00,
     & 0.000000D+00,7.522500D+00,3.109000D-01,0.000000D+00,0.000000D+00/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=41,50) /
C   Nb    Z=41
     & 6.854655D+02,9.699750D+01,1.726960D+01,2.557800D+00,1.950000D-01,
     & 0.000000D+00,8.884600D+01,1.410410D+01,1.577200D+00,0.000000D+00,
     & 0.000000D+00,8.352000D+00,2.677000D-01,0.000000D+00,0.000000D+00,
C   Mo    Z=42
     & 7.212365D+02,1.028863D+02,1.861960D+01,2.793700D+00,2.004000D-01,
     & 0.000000D+00,9.448000D+01,1.532190D+01,1.753600D+00,0.000000D+00,
     & 0.000000D+00,9.319500D+00,3.122000D-01,0.000000D+00,0.000000D+00,
C   Tc    Z=43
     & 7.580685D+02,1.090976D+02,2.015860D+01,3.172600D+00,2.350000D-01,
     & 0.000000D+00,1.004341D+02,1.672630D+01,2.059400D+00,0.000000D+00,
     & 0.000000D+00,1.047110D+01,4.807000D-01,0.000000D+00,0.000000D+00,
C   Ru    Z=44
     & 7.955255D+02,1.151721D+02,2.142740D+01,3.269200D+00,2.094000D-01,
     & 0.000000D+00,1.062531D+02,1.786240D+01,2.113400D+00,0.000000D+00,
     & 0.000000D+00,1.135700D+01,4.011000D-01,0.000000D+00,0.000000D+00,
C   Rh    Z=45
     & 8.340450D+02,1.215701D+02,2.288650D+01,3.510300D+00,2.132000D-01,
     & 0.000000D+00,1.123930D+02,1.918660D+01,2.297700D+00,0.000000D+00,
     & 0.000000D+00,1.242840D+01,4.460000D-01,0.000000D+00,0.000000D+00,
C   Pd    Z=46
     & 8.733135D+02,1.279657D+02,2.420900D+01,3.587300D+00,0.000000D+00,
     & 0.000000D+00,1.185306D+02,2.037420D+01,2.330100D+00,0.000000D+00,
     & 0.000000D+00,1.336340D+01,3.360000D-01,0.000000D+00,0.000000D+00,
C   Ag    Z=47
     & 9.138315D+02,1.348775D+02,2.591770D+01,4.001500D+00,2.200000D-01,
     & 0.000000D+00,1.251811D+02,2.194530D+01,2.676800D+00,0.000000D+00,
     & 0.000000D+00,1.467810D+01,5.374000D-01,0.000000D+00,0.000000D+00,
C   Cd    Z=48
     & 9.553115D+02,1.420059D+02,2.770850D+01,4.450500D+00,2.649000D-01,
     & 0.000000D+00,1.320469D+02,2.359710D+01,3.053500D+00,0.000000D+00,
     & 0.000000D+00,1.607190D+01,7.636000D-01,0.000000D+00,0.000000D+00,
C   In    Z=49
     & 9.977965D+02,1.493945D+02,2.962450D+01,4.976700D+00,3.727000D-01,
     & 0.000000D+00,1.391718D+02,2.537420D+01,3.507200D+00,1.973000D-01,
     & 0.000000D+00,1.758950D+01,1.063100D+00,0.000000D+00,0.000000D+00,
C   Sn    Z=50
     & 1.041224D+03,1.569813D+02,3.160340D+01,5.517000D+00,4.793000D-01,
     & 0.000000D+00,1.464938D+02,2.721350D+01,3.973500D+00,2.486000D-01,
     & 0.000000D+00,1.916790D+01,1.373500D+00,0.000000D+00,0.000000D+00/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=51,60) /
C   Sb    Z=51
     & 1.085597D+03,1.647690D+02,3.364800D+01,6.074800D+00,5.889000D-01,
     & 0.000000D+00,1.540157D+02,2.911790D+01,4.456200D+00,2.992000D-01,
     & 0.000000D+00,2.080990D+01,1.699400D+00,0.000000D+00,0.000000D+00,
C   Te    Z=52
     & 1.130916D+03,1.727580D+02,3.575830D+01,6.650500D+00,7.027000D-01,
     & 0.000000D+00,1.617379D+02,3.108750D+01,4.956000D+00,3.506000D-01,
     & 0.000000D+00,2.251590D+01,2.041800D+00,0.000000D+00,0.000000D+00,
C    I    Z=53
     & 1.177181D+03,1.809481D+02,3.793420D+01,7.244300D+00,8.211000D-01,
     & 0.000000D+00,1.696601D+02,3.312210D+01,5.473300D+00,4.032000D-01,
     & 0.000000D+00,2.428550D+01,2.401200D+00,0.000000D+00,0.000000D+00,
C   Xe    Z=54
     & 1.224392D+03,1.893390D+02,4.017540D+01,7.856200D+00,9.444000D-01,
     & 0.000000D+00,1.777822D+02,3.522140D+01,6.008300D+00,4.573000D-01,
     & 0.000000D+00,2.611870D+01,2.777900D+00,0.000000D+00,0.000000D+00,
C   Cs    Z=55
     & 1.272764D+03,1.981426D+02,4.269280D+01,8.695400D+00,1.231600D+00,
     & 1.237000D-01,1.863159D+02,3.759570D+01,6.768500D+00,6.835000D-01,
     & 0.000000D+00,2.822600D+01,3.379500D+00,0.000000D+00,0.000000D+00,
C   Ba    Z=56
     & 1.322088D+03,2.071532D+02,4.528060D+01,9.556300D+00,1.512700D+00,
     & 1.575000D-01,1.950556D+02,4.003950D+01,7.549300D+00,9.039000D-01,
     & 0.000000D+00,3.040210D+01,4.001500D+00,0.000000D+00,0.000000D+00,
C   La    Z=57
     & 1.372274D+03,2.162757D+02,4.784400D+01,1.034540D+01,1.704400D+00,
     & 1.704000D-01,2.039070D+02,4.245880D+01,8.258900D+00,1.049400D+00,
     & 0.000000D+00,3.255340D+01,4.553600D+00,2.688000D-01,0.000000D+00,
C     Z=58
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=59
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=60
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=61,70) /
C     Z=61
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=62
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=63
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=64
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=65
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=66
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=67
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=68
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=69
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C     Z=70
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00/
      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=71,80) /
C     Z=71
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
     & 0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,0.000000D+00,
C   Hf    Z=72
     & 2.231802D+03,3.693353D+02,8.541780D+01,1.780350D+01,2.524800D+00,
     & 2.104000D-01,3.529815D+02,7.785830D+01,1.462780D+01,1.537100D+00,
     & 0.000000D+00,6.395470D+01,8.886900D+00,2.992000D-01,1.436200D+00,
C   Ta    Z=73
     & 2.296837D+03,3.813025D+02,8.863930D+01,1.869490D+01,2.729700D+00,
     & 2.197000D-01,3.646764D+02,8.092800D+01,1.543500D+01,1.696600D+00,
     & 0.000000D+00,6.674300D+01,9.530900D+00,3.516000D-01,1.815200D+00,
C    W    Z=74
     & 2.362823D+03,3.934733D+02,9.193230D+01,1.961100D+01,2.933900D+00,
     & 2.277000D-01,3.765739D+02,8.406830D+01,1.626580D+01,1.856300D+00,
     & 0.000000D+00,6.960040D+01,1.019650D+01,4.029000D-01,2.214300D+00,
C   Re    Z=75
     & 2.429758D+03,4.058468D+02,9.529570D+01,2.055070D+01,3.138400D+00,
     & 2.347000D-01,3.886733D+02,8.727820D+01,1.711930D+01,2.017100D+00,
     & 0.000000D+00,7.252590D+01,1.088300D+01,4.538000D-01,2.633300D+00,
C   Os    Z=76
     & 2.497642D+03,4.184224D+02,9.872830D+01,2.151300D+01,3.344000D+00,
     & 2.409000D-01,4.009741D+02,9.055660D+01,1.799440D+01,2.179600D+00,
     & 0.000000D+00,7.551840D+01,1.158980D+01,5.048000D-01,3.071700D+00,
C   Ir    Z=77
     & 2.566475D+03,4.311994D+02,1.022294D+02,2.249700D+01,3.550900D+00,
     & 2.466000D-01,4.134756D+02,9.390270D+01,1.889050D+01,2.344000D+00,
     & 0.000000D+00,7.857720D+01,1.231620D+01,5.562000D-01,3.529200D+00,
C   Pt    Z=78
     & 2.636075D+03,4.440131D+02,1.056345D+02,2.333780D+01,3.606700D+00,
     & 2.180000D-01,4.260137D+02,9.715270D+01,1.964400D+01,2.373300D+00,
     & 0.000000D+00,8.153900D+01,1.290050D+01,4.764000D-01,3.843400D+00,
C   Au    Z=79
     & 2.706793D+03,4.571816D+02,1.092606D+02,2.435370D+01,3.808800D+00,
     & 2.208000D-01,4.389052D+02,1.006224D+02,2.057060D+01,2.534300D+00,
     & 0.000000D+00,8.471870D+01,1.365540D+01,5.210000D-01,4.328300D+00,
C   Hg    Z=80
     & 2.778639D+03,4.707342D+02,1.131364D+02,2.557330D+01,4.182000D+00,
     & 2.610000D-01,4.521794D+02,1.043406D+02,2.169890D+01,2.850800D+00,
     & 0.000000D+00,8.814520D+01,1.460960D+01,7.142000D-01,5.012400D+00/

      DATA ((ENERGY(ISHELL,IZ),ISHELL=1,15),IZ=81,86) /
C   Tl    Z=81
     & 2.851505D+03,4.845593D+02,1.171502D+02,2.688340D+01,4.618700D+00,
     & 3.611000D-01,4.657258D+02,1.081966D+02,2.291760D+01,3.231300D+00,
     & 1.924000D-01,9.170840D+01,1.565290D+01,9.683000D-01,5.785200D+00,
C   Pb    Z=82
     & 2.925331D+03,4.985974D+02,1.212423D+02,2.822400D+01,5.059600D+00,
     & 4.589000D-01,4.794846D+02,1.121304D+02,2.416650D+01,3.614700D+00,
     & 2.398000D-01,9.534800D+01,1.672500D+01,1.224500D+00,6.586600D+00,
C   Bi    Z=83
     & 3.000119D+03,5.128515D+02,1.254155D+02,2.959810D+01,5.508200D+00,
     & 5.582000D-01,4.934588D+02,1.161446D+02,2.544830D+01,4.005000D+00,
     & 2.862000D-01,9.906700D+01,1.782880D+01,1.487400D+00,7.419400D+00,
C   Po    Z=84
     & 3.075872D+03,5.273220D+02,1.296699D+02,3.100570D+01,5.965200D+00,
     & 6.602000D-01,5.076490D+02,1.202397D+02,2.676310D+01,4.403200D+00,
     & 3.327000D-01,1.028656D+02,1.896460D+01,1.758300D+00,8.283900D+00,
C   At    Z=85
     & 3.152587D+03,5.420085D+02,1.340053D+02,3.244670D+01,6.431000D+00,
     & 7.654000D-01,5.220545D+02,1.244153D+02,2.811080D+01,4.809800D+00,
     & 3.799000D-01,1.067438D+02,2.013210D+01,2.037900D+00,9.180000D+00,
C   Rn    Z=86
     & 3.230265D+03,5.569115D+02,1.384216D+02,3.392060D+01,6.905800D+00,
     & 8.740000D-01,5.366755D+02,1.286712D+02,2.949110D+01,5.225200D+00,
     & 4.280000D-01,1.107011D+02,2.133120D+01,2.326300D+00,1.010760D+01/
C

      DATA ((ENERGY7(ISHELL,IZ),ISHELL=1,19),IZ=87,103) /
C         Z=87:89
     & 57*0.0D0,
C   Th    Z=90
     & 3.630265D+03,5.969115D+02,1.544216D+02,3.792060D+01,8.905800D+00,
     & 1.274000D+00,2.317762D-01,
     & 5.766755D+02,1.446712D+02,3.349110D+01,6.825200D+00,6.280000D-01,
     & 1.439827D-01,
     & 1.267011D+02,2.533120D+01,3.526300D+00,3.659896D-01,
     & 1.410760D+01,6.131576D-01,
C   Pa    Z=91
     & 3.730265D+03,6.069115D+02,1.584216D+02,3.892060D+01,9.405800D+00,
     & 1.374000D+00,2.164518D-01,
     & 5.866755D+02,1.486712D+02,3.449110D+01,7.225200D+00,6.780000D-01,
     & 1.359714D-01,
     & 1.307011D+02,2.633120D+01,3.826300D+00,2.972997D-01,
     & 1.510760D+01,7.485778D-01,
C    U    Z=92
     & 3.830265D+03,6.169115D+02,1.624216D+02,3.992060D+01,9.905800D+00,
     & 1.474000D+00,2.276235D-01,
     & 5.966755D+02,1.526712D+02,3.549110D+01,7.625200D+00,7.280000D-01,
     & 1.426965D-01,
     & 1.347011D+02,2.733120D+01,4.126300D+00,3.349674D-01,
     & 1.610760D+01,8.589352D-01,
C   Np    Z=93
     & 3.930265D+03,6.269115D+02,1.684216D+02,4.092060D+01,1.040580D+01,
     & 1.574000D+00,2.302694D-01,
     & 6.066755D+02,1.566712D+02,3.649110D+01,8.225200D+00,7.780000D-01,
     & 1.440930D-01,
     & 1.387011D+02,2.833120D+01,4.426300D+00,3.396713D-01,
     & 1.710760D+01,9.281336D-01,
C   Pu    Z=94
     & 4.030265D+03,6.369115D+02,1.724216D+02,4.192060D+01,1.090580D+01,
     & 1.674000D+00,2.214497D-01,
     & 6.166755D+02,1.606712D+02,3.749110D+01,8.625200D+00,8.280000D-01,
     & 1.410060D-01,
     & 1.427011D+02,2.933120D+01,4.726300D+00,2.324009D-01,
     & 1.810760D+01,8.401564D-01,
C         Z=95:103
     & 171*0.0D0/
C
      IF ((NQ.LE.0) .OR. (NQ.GE.104)) GOTO 999
C
C     Calculate offset for this l quantum number for ECP,
C     where the lowest orbitals are excluded.
C     If not ECP, i.e. IQCORE .eq. 0, then ILOFF will be zero.
C
      CALL ECP_LCORE(IQCORE,IQM,ILOFF)
      ISHL = ISHL1 + ILOFF
      IF (NQ .LE. 86) THEN
         HUCTER = -ENERGY(IOFF(IQM)+ISHL,NQ)
      ELSE
         HUCTER = -ENERGY7(IOFF7(IQM)+ISHL,NQ)
      END IF
      IF (HUCTER.EQ.0.0D0) GOTO 999
      RETURN
C
 999  CONTINUE
      WRITE (LUPRI,'(/A,I3/A,4I10//A)')
     &' No orbital ionization potential for nucleus with charge ', NQ,
     &' Debug info: IQCORE, IQM, ISHL, ILOFF =',IQCORE,IQM,ISHL,ILOFF,
     &' INFO: Hueckel start guess replaced with a diagonlization of H1'
      DOHUCKEL = .FALSE.
      RETURN
      END
C   /* Deck hucort */
      SUBROUTINE HUCORT(CMO,CMO1,S,WORK,LWORK)
C
C     Written by K.Ruud, aug-96
C     Version written on the basis of the ORTHO and NORM subroutines
C     of hjj
C     Revised by hjaaj May 2000 to handle round-off errors
C
#include "implicit.h"
#include "maxorb.h"
      PARAMETER (THNORM = 1.0D-6, D1 = 1.0D0, D0 = 0.0D0, THRRND=0.9D0)
      DIMENSION CMO(*), CMO1(*), S(*), WORK(LWORK)
C
#include "priunit.h"
#include "huckel.h"
#include "infinp.h"
#include "inforb.h"
#include "infpri.h"
C
      CALL QENTER('HUCORT')
C
C     Read overlap matrix from LUONEL
C
      CALL RDONEL('OVERLAP ',.TRUE.,S,NNBAST)
      IF (P6FLAG(6)) THEN
         WRITE(LUPRI,'(/A)')
     *     ' (HUCORT) Overlap between AO basis functions :'
         CALL OUTPKB(S,NBAS,NSYM,-1,LUPRI)
      END IF
C
C     We orthogonalize orbitals symmetry by symmetry
C
      ISTS   = 1
      ISTC   = 1
      DO 100 ISYM = 1, NSYM
         NORBI = NORB(ISYM)
         NBASI = NBAS(ISYM)
         NHUCI = NHUCAO(ISYM)
         IF (NORBI .EQ. 0) GO TO 101
         IF (NHUCI .GT. NORBI) THEN
            NOCCI = NOCC(ISYM)
            IF (NOCCI .EQ. 0) NOCCI = NORBI
            NHUCI = MIN( NHUCI, NOCCI )
         END IF
C        ... hjaaj Sep07: take care of calculations with minimum basis set
C                         smaller than Huckel, e.g. S atoms without d-orbitals
C            hjaaj Jun10: only project on the occupied orbitals in this
C                         case (because we have seen a problem when
C                         running an STO-3G calculation with lower symmetry
C                         than the real symmetry)
C
C     We orthogonalize the projected Huckel eigenvectors first
C
         IF (NHUCI .GT. 0) THEN
            CALL NORM(S(ISTS),CMO(ISTC),NBASI,NHUCI,WORK,THNORM,IRETUR)
            IF (IRETUR.NE.0) GO TO 5000
         END IF
C
C     Now we add vectors that span the entire space and orthogonalize 
C     them to the previous vectors
Chj-may2000: modified to use the previous set of vectors in CMO1. The
C     original algorithm here could reinsert CMO vectors deleted by
C     DELMO or otherwise.
C     The initial vectors in CMO1 are generated by H1MO, so the NHUCI
C     first will probably have big overlap with the Huckel generated
C     orbitals. Therefore we begin with no. NHUCI+1 and wrap to 1
C     to end with NHUCI.
C
         IKROW1 = NHUCI
         DO 103 IKROW = NHUCI + 1, NORBI
            IVKC  = ISTC + (IKROW -1)*NBASI
 107        IKROW1 = IKROW1 + 1
            IF (IKROW1 .GT. NORBI) IKROW1 = 1
            IVKC1 = ISTC + (IKROW1-1)*NBASI
            CALL DCOPY(NBASI,CMO1(IVKC1),1,CMO(IVKC),1)
C
C           normalize initial vector C_k and find W_k = S*C_k
C
            ITURN = 0
 108        CONTINUE 
            CALL MPAPV(NBASI,S(ISTS),CMO(IVKC),WORK)
            TNORM = DDOT(NBASI,CMO(IVKC),1,WORK,1)
            IF (TNORM .LT. THNORM) THEN
Chj         ... zero vector on input
               IRETUR = -IKROW
               GO TO 5000
            END IF
Chj may2000: normalize input vector
C        (we ignore round-off errors as it is renormalized later)
            TNORM = D1 / SQRT(TNORM)
            CALL DSCAL(NBASI,TNORM,CMO(IVKC),1)
            CALL MPAPV(NBASI,S(ISTS),CMO(IVKC),WORK)
            TNORM = DDOT(NBASI,CMO(IVKC),1,WORK,1)
C
            IVJC = ISTC
            DO J = 1, IKROW-1
               T = DDOT(NBASI,CMO(IVJC),1,WORK,1)
               TNORM = TNORM - T*T
               IVJC = IVJC + NBASI
               WORK(NBASI+J) = -T
            END DO
            IF (TNORM .LT. THNORM) THEN
C
C              Zero norm,
C              we try the next mo vector from the initial mo's in CMO1
C
               IF (IKROW1 .EQ. NHUCI) THEN
C                 ... we have tried them all!
                  IRETUR = -10 000*IKROW
                  GO TO 5000
               END IF
               GOTO 107
            END IF
C
            IF (TNORM .LT. THRRND) THEN
               IF (ITURN .GT. 1) THEN
                  IRETUR = 10 000*IKROW
                  GO TO 5000
               END IF
               IROUND = 1
            ELSE
               IROUND = 0
            END IF
C
            TNORM = D1/SQRT(TNORM)
            DO J = 1, IKROW-1
               WORK(NBASI+J) = WORK(NBASI+J)*TNORM
            END DO
            WORK(NBASI+IKROW) = TNORM
            CALL DGEMM('N','N',NBASI,1,IKROW,1.D0,
     &                 CMO(ISTC),NBASI,
     &                 WORK(NBASI+1),IKROW,0.D0,
     &                 WORK,NBASI)
            CALL DCOPY(NBASI,WORK,1,CMO(IVKC),1)
C
            IF (IROUND .EQ. 1) THEN
               ITURN = ITURN + 1
               GO TO 108
            END IF
C
 103     CONTINUE 
 101     ISTS   = ISTS   + NBASI*(NBASI + 1)/2
         ISTC   = ISTC   + NORBI*NBASI
 100  CONTINUE 
      CALL QEXIT('HUCORT')
      RETURN
C
 5000 CONTINUE
      WRITE(LUW4,4010) ISYM, IRETUR
      IF (LUPRI.NE.LUW4)
     *   WRITE(LUPRI,4010) ISYM, IRETUR
      WRITE(LUERR,4010) ISYM, IRETUR
      CALL QUIT('*** ERROR *** FATAL ERROR IN HUCORT')
 4010 FORMAT(/' *** HUCORT-FATAL ERROR *** Linear dependency in',
     *        ' symmetry (Huckel part) =',I3/' CODE =',I10)
      END
